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Abstract 

For a projection-based reduced order model (ROM) of a fluid flow to be stable and 
accurate, the dynamics of the truncated subspace must be taken into account. This 
paper proposes an approach for stabilizing and enhancing projection-based fluid ROMs 
in which truncated modes are accounted for a priori via a minimal rotation of the 
projection subspace. Attention is focused on the full non-linear compressible Navier- 
Stokes equations in specific volume form as a step toward a more general formulation 
for problems with generic non-linearities. Unlike traditional approaches, no empirical 
turbulence modeling terms are required, and consistency between the ROM and the full 
order model from which the ROM is derived is maintained. Mathematically, the approach 
is formulated as a trace minimization problem on the Stiefel manifold. The reproductive 
as well as predictive capabilities of the method are evaluated on several compressible flow 
problems, including a problem involving laminar flow over an airfoil with a high angle of 
attack, and a channel-driven cavity flow problem. 

Keywords: Projection-based reduced order model (ROM), Proper Orthogonal 
Decomposition (POD), compressible flow, stabilization, trace minimization, Stiefel 
manifold. 


1. Introduction 

The past several decades have seen an exponential growth of computer processing 
speed and memory capacity. The massive, complex simulations that run on supercom¬ 
puters allow exploration of fields for which physical experiments are too impractical, 
hazardous, and/or costly. A striking number of these fields require computational fluid 
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dynamics (CFD) models and simulations. Accurate and efficient CFD simulations are 
critical to many defense, climate and energy missions, e.g., simulations aimed to de¬ 
sign and qualify nuclear weapons components carried within an aircraft weapons bay; 
global climate simulations aimed to predict anticipated twenty-first century sea-level rise; 
aero-elastic simulations for optimal design of wind systems for power generation. 

Unfortunately, even with the aid of massively parallel next-generation computers, 
CFD simulations for applications such as these are still too expensive for real-time and 
multi-query applications such as uncertainty quantification (UQ), optimization and con¬ 
trol design. Reduced order modeling is a promising tool for bridging the gap between 
high-fidelity, and real-time simulations/UQ. Reduced order models (ROMs) are derived 
from a sequence of high-fidelity simulations but have a much lower computational cost. 
Hence, ROMs can enable real-time simulations of complex systems for more rapid anal¬ 
ysis, control and decision-making in the presence of uncertainty. 

Most existing ROM approaches are based on projection. In projection-based re¬ 
duced order modeling, the state variables are approximated in a low-dimensional sub¬ 
space. There exist a number of approaches for calculating this low-dimensional subspace, 
e.g.. Proper Orthogonal Decomposition (POD) [SI [19], Dynamic Mode Decomposition 
(DMD) [351111] balanced POD (BPOD) [331 [S], balanced truncation [T71[37], and the 
reduced basis method [331117]. In all of these methods, a basis for the low-dimensional 
subspace is obtained from a basis for a higher-dimensional subspace through truncation 
- the removal of modes that are believed to be unimportant in representing a prob¬ 
lem solution. Typically, the size of the reduced basis is chosen according to an energy 
criterion: modes with low energy are discarded, so that the reduced basis subspace is 
spanned by the highest energy modes. Although truncated modes are negligible from 
a data compression point of view, they are often crucial for representing solutions to 
the dynamical flow equations. Dynamics of the truncated orthogonal subspace must be 
taken into account for to ensure stability and accuracy of the ROM. 

For linear systems, a variety of techniques for generating low-dimensional projection- 
based ROMs with rigorous stability guarantees and accuracy bounds are available nzi 
1371 HJ [33]. Equivalent results are lacking for nonlinear systems. Traditionally, low¬ 
dimensional ROMs of fluid flows have been stabilized and enhanced using empirical 
turbulence models. In this approach, the nonlinear dynamics of the truncated subspace 
are modeled using additional constant and linear terms in the ROM system of ordinary 
differential equations (ODEs) [3] [311 [T31 ITl] . More recently, nonlinear eddy-viscosity 
models have also been proposed [3311351 [331133]. One downside of turbulence models 
is that they destroy consistency between the Navier-Stokes partial differential equations 
(PDEs) and the ODE system of the ROM. Accurately identifying and matching free 
coefficients of the turbulence models is another challenge. Moreover, these methods are 
usually limited to the incompressible Navier-Stokes equations. 

Consider, for concreteness, the POD/Galerkin approach to model reduction applied 
to the incompressible Navier-Stokes equations. For these equations, the natural choice 
of inner product for the Galerkin projection step of the model reduction procedure is the 

inner product. This is because, in these models, the solution vector is taken to be the 
velocity vector u, so that ||u ||2 is a measure of the global kinetic energy in the domain, 
and the POD modes optimally represent the kinetic energy present in the ensemble from 
which they were generated. The same is not true for the compressible Navier-Stokes 
equations. Hence, if a compressible fluid ROM is constructed in the inner product (a 
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common choice of inner product in projection-based model reduction), the ROM solution 
may not satisfy the conservation relation implied by the governing equations, and may 
exhibit non-physical instabilities m- 

Unfortunately, ROM instability is a real problem for many compressible flow prob¬ 
lems: as demonstrated in PEiiTiEg, a compressible fluid POD/Galerkin ROM might 
be stable for a given number of modes, but unstable for other choices of basis size. 
Several researchers have proposed ways to circumvent this difficulty through the care¬ 
ful construction of an energy-based inner product for the projection step of the model 
reduction. Rowley et al. m show that Galerkin projection preserves the stability of 
an equilibrium point at the origin if the ROM is constructed in an energy-based in¬ 
ner product. Barone et al. [7], Kalashnikova et al. [22] demonstrate that a symmetry 
transformation leads to a stable formulation for a Galerkin ROM for the linearized com¬ 
pressible Euler equations and nonlinear compressible Navier-Stokes equations with solid 
wall and far-held boundary conditions. Serre et al. [42] propose applying the stabilizing 
projection developed by Barone et al. [7], Kalashnikova et al. [22] to a skew-symmetric 
system constructed by augmenting a given linear system with its adjoint system. The 
downside to these methods is that they are inherently embedded methods: access to the 
governing PDEs and/or the code that discretizes these PDEs is required. 

Other ROM approaches, e.g., the Gauss-Newton with Approximated Tensors (GNAT) 
method of Carlberg et al. m, have better stability properties, as they formulate the 
ROM at the fully discrete level. The drawback of this approach is that an additional 
layer of approximation - usually called hyper-reduction — is required to gain computa¬ 
tional speed-up. Moreover, the approach lacks stability guarantees for low-dimensional 
expansions. 

In this paper, a stabilization and enhancement approach to ROMs for the compress¬ 
ible Navier-Stokes equations is developed. The approach is an extension of the method¬ 
ology developed in mm specifically for the incompressible Navier-Stokes equations. The 
specific volume (C-) form of the compressible Navier-Stokes equations is utilized. Since 
these equations have polynomial (quadratic) nonlinearities, the Galerkin projection can 
be computed offline, once and for all; no hyper-reduction is required m- Unlike tra¬ 
ditional eddy-viscosity-based stabilization methods, the proposed approach requires no 
additional empirical turbulence modeling terms - truncated modes are accounted for a 
priori via a minimal rotation of projection subspace. The method is also non-intrusive, 
as it operates only on the matrices and tensors defining a ROM ODE system, which 
are stabilized through the offline solution of a small trace minimization problem on the 
Stiefel manifold. The proposed new approach can be interpreted as a combination of 
several previously developed ideas. Following lollo et al. [21], we propose to stabilize 
and enhance projection-based ROMs by modifying the projection subspace in order to 
capture more of the low-energy, but high dissipative scales of the flow solution. Similarly 
to Amsallem and Farhat [T], a rotation of the projection subspace is used to achieve this 
goal. Specifically, a larger set of basis is linearly superimposed to provide a smaller set of 
basis that generate a stable and accurate ROM. Finally, in the spirit of most previously 
proposed eddy-viscosity-based turbulence models, the eigenvalues of the linear part of 
the Galerkin system of ODEs are used as a proxy to guide the stabilization algorithm. 

The remainder of this paper is organized as follows. In § [^ the standard projection- 
based model reduction approach is outlined in the context of the specific volume form 
of the compressible Navier-Stokes equations. Also overviewed is the POD method for 
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constructing an optimal reduced basis, and some eddy-viscosity based closure models 
used for accounting for modal truncation. The proposed methodology of “rotating” the 
projection subspace into a more dissipative regime to better resolve the small, energy 
dissipative scales of the flow is detailed in § [^ Here, the approach is formulated math¬ 
ematically as a constrained optimization problem on the Stiefel manifold. In § [^ the 
performance of the proposed method is evaluated on several compressible flow problems, 
including a problem involving a laminar flow over an airfoil with a high angle of attack, 
and a channel-driven cavity flow problem. Finally, conclusions are offered in § [^ 


2. Projection-based model reduction for nonlinear compressible flow 

In this section the standard projection-based model reduction approach is laid out. 


In § |2.1| the approach for a general nonlinear system is presented while in § 2^ the 
approach is applied to the compressible Navier-Stokes equations. The POD method for 
calculating a reduced basis using a set of snapshots from a high-fidelity simulation is 
outlined in § |2.3[ followed by a brief overview of eddy-viscosity-based closure models 
that account for modes truncated in the application of the POD method (§ 2.4). 


2.1. Nonlinear projection-based model order reduction 
Consider a dynamical system of the form: 

^w = F{w), (1) 

where F is the propagator in H, a Hilbert space. In fluid flows, the state variable 
w — w{x, t) G H depends on space a; S O, D being the flow domain, and time t G [0, T], 
T representing the period of integration. Then, the propagator F contains spatial deriva¬ 
tives. The associated Hilbert space of square-integrable functions L^(0) is equipped with 
the standard inner product for its elements w G defined by: 


{w,w)q:= j w ■ w dx. ( 2 ) 

n 

In the Galerkin ROM approach, the governing variable, w(x^t) is discretized using 
basis functions (modes) ^ H with corresponding mode coefficients 

n 

w{x,t) « Wo(x) -I- w^^--"^w,t) := Wo(x) -t ^ Oj(t) Wi{x), (3) 

z=l 

where Wq{x) denotes the (steady) mean flow. 

In the method of lines, the modes Wi are known a priori and the goal is to find 
mode coefficients Oi that satisfy the differential equation Q. In general, the modes 
Wi can be chosen in a number of ways. In the context of spectral methods in CFD for 
example, the basis vectors are usually analytical functions, e.g. trigonometric functions or 
Chebyshev polynomials. The advantage of these functions is that their spatial derivatives 
have analytical representations and numerically efficient algorithms such as the Fast 
Fourier Transform (FFT) can be utilized. In the context of ROMs, the spatial basis 
functions are usually derived a posteriori from a snapshot of a solution data set, like 
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the Proper Orthogonal Decomposition (POD) [TS] or Dynamic Mode Decomposition 
(DMD) [5S1I3T]. Attention is restricted here to modes computed using the POD method 
(detailed in § 2.3), but it is noted that the methods proposed here hold for any choice 
of reduced basis. The reason for the choice of the POD reduced basis is two-fold. First, 
the POD is a widely used approach for computing efficient bases for fluid dynamical 
systems. Moreover, ROMs constructed via the POD/Galerkin method lack in general an 
a priori stability guarantee (meaning POD/Galerkin ROMs would benefit from ROM 
stabilization approaches such as the one developed herein). 

The mode coefficients in ^ Oi are chosen to minimize the residual of the Galerkin 
expansion 


for i = 1,..., n. This projection yields a set of evolution equations for the mode coefficients 

tti 


d 




= Ma), 


(5) 


where a := (ai,...,a„)^ represents the state and / := (/i,...,/„)^ its propagator. 
Given some initial conditions, the evolution equation ([^ can be integrated using standard 
numerical integration techniques. The ROM system ([^ is, by construction, small, and 
can be integrated forward in time in real or near-real time unlike the full order dynamical 
system from which it is derived. 


2.2. Nonlinear reduction of the compressible Navier-Stokes equations 

Gonsider the 2D compressible Navier-Stokes equations in primitive variable^ 


Ct T ^Cx T ^Cy 

Ut + UUx + VUy + (Px = :^C 

Re 


Vt + UVx + VVy -I- CPy = 

Ke 


Pt + UPx + VPy + JpiVx +Uy) = 


7 


RePr 
1 - 7 
Re 


^Ux gUj 


^Vy 


[{pOxx + (K) 


+ {Vx -f Uy)y 

T i'^x “t” '^y)x 


( 6 a) 

( 6 b) 

( 6 c) 


Ux ( yMx - y^’^/ ) + I’y ( ^Vy - -Ux ) + (Wy -|- VxY 

( 6 d) 


yyi 


Here, ({x, t) = Xj p{x, t) is the specihc volume (the inverse of the density, p{x, t)), u{x, t) 
and v{x,t) are the Gartesian components of the flow velocity, p{x,t) is the pressure, 7 
is the specific heat ratio. Re is the Reynolds number, Pr is the Prandtl number, and the 


^Presented in two-dimensions for the sake of brevity only. Extension to the three-dimensional equa- 
tions is straightforward. 
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subscripts denote partial derivatives. A Galerkin projection yields a system of coupled 
quadratic ODEs whose constant coefficients are calculated off-line and once and for all 
and lollo et al. [H] for details). This system has the form: 

— C + La 'y a ■■■ , (7) 

where C G K", i G and G Vi = 1,..., n. 

Remark 1: The compressible Navier-Stokes equations are typically expressed in conser¬ 
vative form. This form is convenient for many applications including CFD. The conser¬ 
vative form contains rational functions of the unknowns and it is therefore not possible 
to pre-compute ROMs using standard Galerkin projection; to attain any computational 
speed-up a hyper-reduction step is necessary. Hyper-reduction is not always desirable, 
as it can destroy energy conservation properties and/or symplectic time-evolution maps 
nnna. On the other hand, if the equations are expressed in primitive variables, hyper¬ 
reduction can be avoided because all nonlinearities that appear are polynomial. It is for 
this reason that, in our approach, we base the ROM on the equations (§. It is possible 
to extend our proposed approach to problems where the full-order model (FOM) is based 
on the compressible Navier-Stokes equations in conservative form; see Remark 2. 

2.3. Construction of optimal reduced-order basis via the POD 

As discussed earlier, there exist a number of methods for calculating a reduced basis 
G H, e.g., proper orthogonal decomposition (POD) [3H|TH], Dynamic Mode 
Decomposition (DMD) [38l|4T] balanced POD (BPOD) [MlllQ], balanced truncation [13 
nn, and the reduced basis method [33 EH- In this paper, attention is restricted to reduced 
bases constructed using the first of these approaches, namely the POD method. This 
method is reviewed succinctly below. 

Discussed in detail in Lumley [5S] and Holmes et al. [13, POD is a mathematical 
procedure that, given an ensemble of data and an inner product, constructs a basis for 
the ensemble. The POD basis is optimal in the sense that it describes more energy (on 
average) of the ensemble in the chosen inner product than any other linear basis of the 
same dimension n. Let to" G denote a snapshot vector, computed as the solution 
of the fully discretized version of Eq. ([^. Here, N is the number of grid points in the 
full order model, for some instance of its parameters — that is, for some specific time 
t, some specific value of the set of flow parameters, or some boundary/initial conditions 
underlying this governing equation. Suppose a total of AT G N snapshots are collected 
from a high fidelity simulation. A snapshot matrix is defined as a matrix M G 
whose columns are individual snapshots. The main focus of this paper is on unsteady 
flows and on snapshots associated with different time-instances. Hence, M.^i := for 
i = 1,...,K. Mathematically, POD seeks an n-dimensional (n ^ K and n ^ 4A^) 
subspace spanned by the set such that the difference between the ensemble 

and its projection onto the reduced subspace is minimized on average. That is, 
a POD basis is obtained by solving the following low-rank matrix approximation problem: 

For a given snapshot matrix M G , find a lower rank matrix M G 

solves the minimization problem 


(see Appendix A 


da 

dt 
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(8) 


where n <C 47V. 


min ||Af —Afll 


rank(A<f)=n 




In this problem, the rank constraint can be taken care of by representing the unknown 
matrix as M = UV, where U S ]^47Vxn y ^ problem ^ becomes 

min \\M-UV\\f. (9) 

C/GR<‘”x",VGK"X-ff 

It is well-known that the solution of the above low-rank approximation problem is given 
by the Eckart-Young-Mirsky [131 US] theorem via the Singular Value Decomposition 
(SVD) of M. Specifically, U = U,^i,n and V = ^here X = U'EV'^. This is 

the so-called “method of snapshots” for computing a POD basis [44] . 

2.4- Aceounting for modal truncation: eddy-viscosity based closure models 

In the Kolmogorov description of the turbulence cascade, the large, energy-carrying 
flow scales transfer energy to successively smaller scales where finally dissipative forces 
can dissipate their energy [ 321131 ] ■ The large, energy-carrying scales are associated with 
the large singular values of the snapshot matrix M, while the smaller, energy-dissipative 
scales of the flow are associated with the smaller singular values. Since low order POD- 
based ROMs remove modes corresponding to small singular values, these ROMs are, by 
construction, not endowed with the dissipative dynamics of the flow. 

Many of the popular methods for accounting for truncated modes fall in to the family 
of eddy-viscosity based closure models. In this family of methods, dynamics of the 
truncated modes are modeled by modifying the constant and linear terms of the Galerkin 
model. In the most general case, the resulting modified Galerkin system is of the following 
form: 


-—= CLa\ ••• F . ( 10 ) 

dt ■' 

For a detailed review of the performance of the various methods based on this ap¬ 
proach the reader is referred to Wang et al. |48| . Broadly speaking, the modified constant 
and linear terms, C, and X, are constructed in an effort to decrease energy production and 
increase energy dissipation. For example, consider the “modal constant eddy-viscosity” 
approach first proposed by Rempfer and Fasel [Mj. In this approach, the linear term is 
modified using a constant matrix, L — L L and the constant term is unmodified. This 
approach may be interpreted as a modification of the eigenvalues of the linear opera¬ 
tor. In other words, the additional linear term is constructed to decrease the magnitude 
of the real positive eigenvalues of L (i.e., decrease energy production) and increase the 
magnitude of the negative eigenvalues of L (i.e., increase energy dissipation)]^ Since 


®Other ROM stabilization approaches, developed independently from the “modal constant eddy- 
viscosity” approach and for a broader range of applications than fluid mechanics, e.g., the method of 
Kalashnikova et al. m for stabilizing generic linear ROMs via optimization-based eigenvalue reassign¬ 
ment, give rise to a similar modification to L. 
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the appropriate eigenvalue distribution is not known a priori, the additional linear term 
must be identified by solution matching techniques. 

The linear eddy-viscosity approach has been successfully applied to a large number of 
Galerkin models of complex, high-Reynolds number flows. The method has been partic¬ 
ularly successful for ROMs of the incompressible Navier-Stokes equations where energy 
conservation and symmetry arguments can be used to further motivate the approach. 

The main drawback of the eddy-viscosity based stabilization approach is the loss of 
consistency between the Navier-Stokes equations and the Galerkin system. The Galerkin 
system is modified empirically and, therefore, the quadratic system of ODEs no longer 
corresponds to a Galerkin projection of the Navier-Stokes equations. In the following 
section, a novel stabilization and enhancement approach that retains consistency is in¬ 
troduced. 

3. Stabilization and enhancement of compressible flow ROMs via subspace 
rotation 

In the new approach proposed herein, the truncated modes are modeled a priori by 
“rotating” the projection subspace into a more dissipative regime. Specifically, instead of 
approximating the solution using only the first n most energetic POD modes, the solution 
is approximated using a linear-superposition of n -I- p (with p > 0) most energetic POD 
modes. Mathematically this can be expressed as: 

n+p 

Wi = '^XjiWj i = l,---,n, (11) 

1=1 

where X G is the orthonormal (X^X = Inxn) “rotation” matrix. The Galerkin 

system tensors associated with these new modes are expressed as a function of X as 
follows: 


n+p 

Qfk = E i, 3 ,k=l,--- ,n, (12a) 

s,q,r—l 

L = X'^LX, (12b) 

C = X'^C*, (12c) 

where C G R”+-P, L G K(”+p)^("+p) and G IR("+p)x("+p)^v 1 = 1, • • • , (n -|-p) are the 
Galerkin system coefficients corresponding to the first n-|-p most energetic POD modes. 
The new Galerkin system is hence: 

— = C -I- La + [ ] , (13) 

where the matrices Q^’'\ L and C are given by Q. 

The goal of the proposed approach is to find X such that 1.) the new modes Wi 
remain good approximations of the flow, and 2.) the new Galerkin ROM is stable and 
accurate. To ensure that these properties are satisfied, a constrained optimization prob¬ 
lem is formulated for X. To guarantee that the new modes remain good approximations 
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of the flow, the distance ||-Xi — In+p,n\\F is minimized, where In+p,n are the first n 
columns of an n +p identity matrix, and || • Hi? is the Frobenious norm To ensure that 
the ROM is stable and accurate, an eddy-viscosity-based constraint equation is used as 
a proxy. In other words, the eigenvalues of the linear operator are modified until a stable 
and accurate ROM is generated. 

Mathematically, the constrained optimization problem for X outlined above reads as 
follows: 


where ry € K and 


minimize 

XGV(„+^).„ 

subject to 


fr -^(n-t-p)xn) 

tviX^LX) = r] 


(14) 


'^{n+p),n ^ ^ 


T+p) X r 


= /„ , p > 0}. 


(15) 


In (151, V(„_|_p)_„ is the Stiefel manifold, defined as the set of (n+p) x n matrices satisfying 
the orthonormality condition X"'"X = [321115]. In Equation (1^ the objective function 

is simplified by utilizing the property that for a real matrix | X|||. = tr(A"'’A). Thus, 
minimizing ||X — In+p,n\\F is equivalent to minimizing —tr(X"^/(„_|_p)xn)- 

The appropriate eigenvalue distribution p, must be identified using a solution match¬ 
ing procedure. Discussion of an approach for selecting rj is deferred until § |3.2[ 


Remark 2: In this paper, we assume that the ROM advanced forward in time during the 
online time-integration step of the model reduction is a system of the form (101, which 
arises when projecting the compressible Navier-Stokes equations in primitive specific 
volume form (|^ onto the reduced basis modes. As suggested in Remark 1, the method 
described here can be applied in the case the ROM is based on the compressible Navier- 
Stokes equations in conservative form (with or without hyper-reduction). In this case, 
the model reduction would proceed as follows: 


Step 1: Run a high-fidelity code to generate snapshots from which the POD basis will be 
constructed. 

Step 2: Construct from the snapshots collected in Step 1 a POD basis for the 

primitive variables. 


Step 3: Project the compressible Navier-Stokes equations in primitive specific-volume form 
(§ onto the modes from Step 2 to obtain a system of the form ([To| . 


Step 4- Use the matrices L and to obtain from the original basis a stabi- 

■rn-\-p 


lized basis Wi = 


J % 5 ^ - f 5 


, n, where X is the solution to (14|. 


Step 5: Transform the stabilized basis {wi(x)'\^^^ into conservative variables, and use it in 
a ROM code that projects the compressible Navier-Stokes equations in conservative 
form (with or without hyper-reduction). 


To apply this procedure, two ROM codes are required: a ROM code that projects the 
compressible Navier-Stokes equations in primitive specific-volume form, and a ROM code 
that projects the compressible Navier-Stokes equations in conservative form. The former 
code is only needed to calculate the L and matrices, which are used for the basis 
stabilization. 
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3.1. Solution of constrained optimization problem 

A common method for solving constrained optimization problems of the form (14) 
is the method of Lagrange multipliers [29]. In this method, the Lagrangian of the opti¬ 
mization problem is computed, and its stationary points are sought, yielding necessary 
optimality conditions for local maxima and minima. The reader can verify that the 


Lagrangian for Eq. (14) is 


C{X,A„A2) := -tr(XT/(„+p)^„)+tr(Ai(XTx*X-S/„))+tr(A2(XTX-J„)), (16) 

where Ai and A 2 are diagonal matrices of Lagrange multipliers. 

Suppose that Ai is a local maximizer of problem (14). Then X satisfies the first-order 
optimality condition Cx = —I(n+p)xn + ^liL* + L* )X + 2 A 2 A = 0, iY{X^L*X — 
^In) = 0, and X'^X — = 0. Solving Eq.(|l4| using Lagrange multipliers is possible; 


however, it is inefficient. Signihcant speed-ups are possible by satisfying the orthonormal¬ 
ity constraint directly via optimization on the Stiefel matrix manifold. In this method, 
with the help of the augmented Lagrange method, the constrained optimization problem 
is reduced to an unconstrained optimization problem on the Stiefel manifold as follows: 


minimize - tr(ATj(„+p)x„) + f^tT{X'^L*X - _ XctiiX^L*X - 2/„), 

XgV(„+p),„ / 

(17) 

where fik is increased until the constraint is satisfied to some desired precision. The 
variable Xc is an estimate of the Lagrange multiplier and is updated according to the 
rule 


A£.+i ^ A£, - Mfctr(AW (18) 

where is the solution of the unconstrained problem at the step. In this work, 

the McLuopt AlATLAB toolbox [^ is used to solve All derivatives in the optimization 
algorithm are calculated analytically. 

3.2. Solution matching procedure for identification of stabilizing eigenvalue spectrum rj 
In this work, the distribution of stabilizing eigenvalues 77 in 0 is identified using 
a brute force approach. A brute force approach is feasible in this case because the 
associated ROMs are small and relatively inexpensive to integrate numerically. 

Given the initial guess the optimization problem is solved and the correspond¬ 
ing transformation matrix X^^^ identified. Using Eq. (12 1 the new Galerkin ROM is 
constructed and integrated numerically. A global ROM “energy” is dehned as follows: 
E{t)^^'> = ^ linear fit of the temporal data is performed: E{t)^^^ « 

Cg L If the absolute value of c\ ' is below some predefined tolerance, the ROM is deemed 
stable and the procedure halts. If the value of c\ ' is positive (indicating increasing en¬ 
ergy) then + e, where e < 0. If the value of is negative (indicating 

decreasing energy) then r/O+i) = _|_ g where e > 0. This search is automated using 

MATLAB’s fzero root finding algorithm. The overall stabilization procedure is sum¬ 
marized in Algorithm In Remark 3 below, we give some guidelines for selecting the 
parameters in Algorithm [l] and some general remarks about the numerical solution of 
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the optimization problem (14). 


Remark 3. Here, we provide some general guidelines and remarks pertaining to Algo¬ 
rithm 1 . 


• We recommend initializing the algorithm with = J„ and 77 ^°^ = = 

tr{Li.n,i-.n), which corresponds to the trivial solution. 

• Our numerical experiments suggests n = p provides best performance; however the 
optimal choice of p remains an open question. 


• The proposed algorithm requires that p > 1. For p = 0, the rotation matrix X is 
square and therefore corresponds to a coordinate transformation that, by definition, 
can have no effect on the dynamics of the system. 


Existence and uniqueness results for the solution to (141 are not provided in this 
work. We have found the results to be insensitive to the initial conditions, and the 


solutions to (14) to be unique for a large number of random initial conditions. 


For some choices of p and r] the optimization problem 0 has no solutions. Let 
Ai > A 2 > ■ ■ • > A„+p be the eigenvalues of the symmetric part of L and let V be a 
matrix of the associated eigenvector^ By definition tr(X) = Therefore, 

(Aj) < < Sr=n+i('^i) ^ necessary condition for the existence of a solution 

to (17). The lower and upper bounds for this ry are provided by the transformations 


X = Vj:n and X = V^n+i-.n+p, respectively. 


4. Numerical experiments 

In this section, we evaluate the performance of Algorithm [l] for stabilizing compress¬ 
ible flow ROMs on several 2D problems: a problem involving a laminar flow around an 
inclined airfoil, and a channel driven cavity problem at two Reynolds numbers. In all 
cases, the flow is governed by the full compressible Navier-Stokes equations with constant 
viscosity. Direct Numerical Simulations (DNS) are performed and POD basis functions 
are derived from snapshots collected during these simulations. ROMs are derived by 
projecting the fully compressible Navier-Stokes equations in specific volume (C-) form 
onto the hrst n most energetic basis modes. The projection is performed off-line, and 
once and for all, resulting in a system of n coupled quadratic ODEs in the form of ([^. 
From this point forward, such ROMs are referred to as “standard POD-Galerkin ROMs”, 
where “standard” refers to the fact that no model is used to account for the dynamics of 
the truncated modes, n-|-1, n-|-2, • • • , 00 . ROMs derived using the stabilization approach 
proposed in this paper are referred to as “stabilized POD-Galerkin ROMs”. 


®We can look at the symmetric part with out loss of generality since tr(i) = tr((i + L'^)l2) 
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Algorithm 1: Stabilization algorithm for compressible Navier-Stokes equations 

input ; Initial guess for stabilizing eigenvalue distribution ROM size n and 
p> 1, Galerkin system matrices associated with the first n + p most 
energetic POD modes. C G L G M("+p)^("+p) and 

Q(d g j = 1, • ■ ■ +Pj convergence tolerance TOL, and 

maximum number of iterations, kmax 
output: Stabilizing rotation matrix X 


1 for fc = 0, • • • , k^ax do 

2 Solve constrained optimization problem on Stiefel manifold: 

minimize 
subject to 


-tr 


8 

9 

10 

11 


Construct new Galerkin matrices using (12) 

Integrate numerically new Galerkin system 

Calculate modal energy 

Perform linear fit of temporal data ^ 0 ^^ 

(k) 

Based on energy growth c\ , calculate e using root-finding algorithm and 
perform update -(- e 


-1 -I- Cg ^ 


if 




< TOL then 


X ■= aw 

terminate the algorithm 


end 


12 end 
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4..1. High angle of attack laminar airfoil 

The first test case involves the 2D flow around an inclined NACA0012 airfoil at Mach 
0.7, and Re = 500 at 25 degrees angle of attack. The Reynolds number is based on the 
chord of the airfoil, c. At this Reynolds number and angle of attack, the flow is separated 
and the solution corresponds to a stable limit cycle. High-fidelity simulation snapshots 
are generated using a second-order, finite-difference, embedded boundary (EB) solver. 
The no-slip adiabatic boundary conditions are satisfied using a first-order ghost fluid 
method. For a detailed description of the scheme the reader is referred to Balajewicz 
and Farhat [3]. The snapshots correspond to a DNS of the compressible Navier-Stokes 
equations. The viscosity of the fluid inside the domain is assumed constant. The com¬ 
putational domain extends 20c in all directions and a sponge zone of thickness 2c is used 
to absorb all waves entering and leaving the domain. The domain is discretized using a 
non-uniform 300 x 300 tensor grid. The flow is initialized by setting the solution at all 
grid points to the free stream values. Time integration is performed using an implicit 
second order BDF scheme and a constant time step corresponding to a CFL=1 is used. 
Snapshot collection begins after 5000 time steps to ensure the solution has reached the 
limit cycle. A total of AT = 500 snapshots are collected every 5 simulation time steps. 
The first four basis functions capture approximately 86% of the snapshot energy. 

Figur e[T]illu strates the performance of a stabilized n, p = 4 ROM of the laminar airfoil. 
In Figure 1(a) the global energy of the DNS, standard (i.e., unstabilized) and stabilized 
ROMs are illustrated. The stabilized ROM is shown to track very accurately the mean 
of the fluctuating energy of the original DNS snapshots while the standard ROM over 
predicts the mean by an order of magnitude. To investigate the long term stability of the 
stabilized ROM, the system was numerically integrated 100x the duration of the original 
snapshots. No change or drift in trajectory was observed during this long integration 
period. In Figure [T(b)| the trajectories of the first and second temporal coefficient, ai(t), 
and 02(1) respectively, are illustrated. The stabilized ROM accurately reproduces the 
closed orbit of the stable limit cycle while the standard ROM predicts an unstable spiral. 
Figures 1(a) and |l(b)| demonstrate how the stabilized ROM reproduces reliably both 
the global mean and fluctuating components of the DNS snapshots. The stabilizing 
transformation matrix X for this problem is illustrated Figure |l(c)[ As expected the 
rotation of the projection subspace is small as demonstrated by the fact that X « 
I(n+p)xn- For this configuration, the normalized error defined as — I^n+pix-nWr/n is 
0.083. 

Finally, the predicted velocity magnitude at the final snapshot is illustrated in Fig¬ 
ure]^ The stabilized ROM (Figure [^a)) reproduces the velocity contours of the original 
high-fidelity DNS simulation (Figure [^c)) remarkably well, in contrast to the standard 
ROM (Figurej^b)). This demonstrates the effectiveness of the proposed model reduction 
approach. 


4-2. Channel driven laminar cavity 

For the results presented in this section, the high-fidelity fluid simulation data are 
generated using a Sandia National Laboratories’ in-house finite volume flow solver known 
as SIGMA CFD. This code is derived from LESLIE3D [ID], a Large Eddy Simulations 
(LES) flow solver originally developed in the Computational Combustion Laboratory at 
the Georgia Institute of Technology. The code has LES as well as DNS capabilities. 
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Figure 1: Nonlinear model reduction of the laminar airfoil. Evolution of modal energy (a), and phase 
plot of the first and second temporal basis, ai(t) and a 2 {t) (b); DNS (thick gray line), standard n = 4 
ROM (dashed blue line), stabilized n,p = 4 ROM (solid black line). Stabilizing rotation matrix, X (c) 


For the channel driven laminar cavity problem considered here, the code was run in 
DNS mode. For a detailed description of the schemes and models implemented within 
LESLIE3D, the reader is referred to [T51ITB]. 

ROMs for the channel driven laminar cavity problem are constructed using a Sandia 
in-house parallel model reduction code known as Spirit, which constructs ROMs 

for compressible flow problems using the POD and continuous projection method. This 
code, detailed in |22] , reads in the snapshot and mesh data written by a high-fidelity flow 
solver, creates a finite element representation of the snapshots and computes the numer¬ 
ical quadrature necessary for evaluation of the inner products arising in the Galerkin 
projection step of the model reduction. All calculations are performed in parallel using 
distributed matrix and vector data structures and parallel eigensolvers from the Trilinos 
project [18], which allows for large data sets and a relatively large number of POD modes. 
The libmesh finite element library |24] is used to compute the element quadratures. 

In the discussion that follows, two variants of the 2D channel driven laminar cavity 
problem are considered: a low Reynolds number variant (Re « 1500) and a moderate 
Reynolds number variant (Re « 5500). Both tests cases involve a Mach 0.6 viscous 
laminar flow over a cavity in a T-shaped domain (Figure]^. The flow conditions for 
both tests are similar to case L2 in [3S]. The free stream pressure is 25 Pa, the free 
stream temperature is 300 K, and the free stream velocity is 208.8 m/s. The viscosity p, 
is spatially constant and calculated such that the above Reynolds numbers are achieved. 
The thermal conductivity k is also constant, calculated such that the Prandtl number 
is Pr= 0.72. At the inflow boundary, a value of the velocity and temperature that 
is above the free stream values is specified. The flow at the cavity walls is assumed 
to be adiabatic and to satisfy a no-slip condition. The remaining outflow boundaries 
are open, and a far-held boundary condition that suppresses the rehection of waves 
into the computational domain is implemented here. The high-hdelity simulation is 
initialized by setting the flow in the cavity to have a zero velocity, free stream pressure, 
and temperature. The region above the cavity is initialized to free stream conditions and 
the how is allowed to evolve. As both SIGMA GFD and Spirit are 3D codes, a 2D mesh 
of the domain D is converted to a 3D mesh by extruding the 2D mesh in the t-direction 
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Figure 2: Snapshot of high angle of attack airfoil at final snapshot; contours of velocity magnitude. 
DNS (top), standard n = 4 ROM (middle), and stabilized n,p = 4 ROM (bottom) 
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by one element. Finally, it is noted that SIGMA CFD and Spirit work with different 
meshes. The former code requires a structured hexahedral discretization, whereas the 
latter assumes a tetrahedral discretization. To overcome this difference, the hexahedral 
high-fidelity meshes associated with the snapshots are converted to tetrahedral meshes 
prior to constructing the ROMs. This is accomplished by breaking up each hexahedral 
element into six tetrahedral elements. 


4 ..2.1. Low Reynolds number (Re « 1500) 

For the first, low Reynolds number variant of the channel driven laminar cavity 
problem, the free-stream viscosity is set to ^00 = 3.17 x 10“® kg/(m-s), so that Re = 
1453.9. The discretized domain, illustrated in Figure consists of 98,408 nodes, cast 
as 292,500 tetrahedral finite elements within the ROM code. Spirit. The 2D extent of 
the domain is: [(-6.42,10) x (-1,10)]\[(-6.42,10) x (-1, 0) U (2,10) x (-1,0)] m. The 
reader can observe that the mesh is structured but non-uniform. 

The high-fidelity solver, SIGMA CFD, is initiated with the conditions described above 
and allowed to run until a statistically stationary flow regime is reached. At this point, 
a total of Kmax = 500 snapshots are collected from SIGMA CFD, taken every At snap = 
1 X 10“^ seconds. The snapshots are used to construct a POD basis of size 4 modes in 
the inner product. This basis captured about 91% of the snapshot energy. For more 
details on this test case, the reader is referred to [55]. 

Figure [^illustrates the performance of a stabilized n,p = A ROM of the low Reynolds 
number cavity. In Figure [4(b)[ the modal energy of the DNS, standard, and stabilized 
ROMs are illustrated. The stabilized ROM is shown to track very accurately the en¬ 
ergy of the original DNS snapshots while the standard ROM is unable to reproduce this 
trajectory. The long term stability of the stabilized ROM was validated by numerically 
integrating the system 100 x the duration of the original snapshots. No change or drift 
in trajectory was observed during this long integration period. In Figure [4(b)| the tra¬ 
jectories of the first and second temporal coefhcient, ai(t), and 02 (t) respectively, are 
illustrated. The stabilized ROM predicts correctly the closed orbit of the stable limit 
cycle while the standard ROM predicts an unstable spiral. The stabilizing transforma- 

As before, the rotation of the 


tion matrix X for this problem is illustrated Figure 4(c) 


For this 


projection subspace is small as demonstrated by the fact that X « I(n+p)xr 
configuration, the normalized error defined as ||-X' — I(n+p)xn\\F/n is 0.1182. 

Figure]^ shows the Power Spectral Density (PSD) of the predicted pressure fluctua¬ 
tions at the bottom right corner of the cavity, x = (2,-1), Both the fundamental and 
first harmonic of the response is accurately predicted by the stabilized n,p = A ROM. 
The PSD of the GFD signal was computed using all available snapshots from t = 0 
to t = 380 where t is non-dimensional. On the other hand, the PSD of the stabilized 
ROM was computed from the signal 100 x past the duration of the original snapshots; 
i.e. t = (38000 - 380) to t = 38000. 

Finally, a snapshot of the predicted velocity magnitude at the final snapshot is illus¬ 
trated in Figure]^ The stabilized ROM (Figure |^a)) reproduces the velocity contours 
of the original high-fidelity DNS simulation (Figure [^c)) remarkably well. In contrast, 
the standard ROM (Figure |^c)) is unstable and inaccurate. 
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J^.2.2. Moderate Reynolds number (Re « 5500) 

The next test case considered is also a channel drive laminar cavity problem, but at a 
higher Reynolds number. The only parameter that is different is the free-stream viscosity, 
now set to /ioo = 8.46 x 10“^ kg/(m-s), so that Re = 5452.1. Also changed is the size 
of the geometry extent, which has a larger sponge region near the outflow regions. This 
is needed to suppress adequately the reflection of waves into the computational domain 
for this problem. Toward this effect, the 2D extent of the domain is: [(—6.42,30) x 
(—1, 30)]\[(—6.42, 30) X (—1,0) U (2,30) x (—1,0)] m. The geometry is discretized by 
117,328 nodes, cast as 345,900 tetrahedral elements in Spirit. As before, the mesh is 
structured but non-uniform. The flow is significantly more chaotic than the Re « 1500 
case considered in Section 11.2.11 

A total of AT = 500 snapshots are collected from SIGMA CFD at increments Atgnap = 
1 X 10“® seconds. As before, snapshots collection does not begin until a statistically 
stationary flow regime has been reached. From these snapshots, a POD basis of size 20 
modes is constructed in the inner product. This basis captures about 72% of the 
snapshot energy. Typically, n would be selected such that the POD basis captures a 
greater percentage of the snapshot ensemble energy (e.g., « 90% or more). We choose 
a basis that captures less energy of the snapshot set to highlight the effectiveness of our 
approach for low-dimensional POD expansions. 

Figure illustrates the performance of a stabilized n,p = 20 ROM of the higher 
Reynolds number cavity problem. In Figure [7(a)| the modal energy of the DNS, standard, 
and stabilized ROMs are illustrated. The standard ROM is shown to over predict the 
energy of the original DNS snapshots by an order of magnitude. The predictive power 
of the stabilized ROM is demonstrated by numerically integrating the ROM 10 x the 
duration of the original snapshots. The stabilizing transformation matrix X for this 
problem is illustrated Figure 7(b) As before, the rotation of the projection subspace 


is small as demonstrated by the fact that X « I(^n+p)xn- For this configuration, the 
normalized error defined as jjX - I{n+p)xn\\F/n is 0.0384. 

Figures and shows the PSDs of the predicted pressure fluctuations at locations 
Xi = (2, —0.5) and X 2 = (2, 0.5), respectively. The PSD of the CFD signal was computed 
using all available snapshots from t = 0 to t = 67 where t is non-dimensional. On the 
other hand, the PSD of the stabilized ROM was computed from the signal 10x past 
the duration of the original snapshots; i.e. t = (670 — 67) to t = 670. The stabilized 
ROM accurately predicts the chaotic pressure fluctuations at both locations. Figure [l0| 
illustrates the Cross Power Spectral Density (CPSD) for pressure fluctuations at Xi and 
X 2 ■ Both the power and phase lag at the fundamental frequency, and the first two super 
harmonics (normalized frequency (xtt rad/sample) « 0.18,0.35, and 0.53) are predicted 
accurately using the stabilized ROM. The phase lag at these three frequencies in Figure[T0| 
as predicted by the CFD and the stabilized ROM is identified by red squares and blue 
triangles, respectively. As expected, the low-dimensional ROM is unable to reproduce 
the phase lag of low-amplitude frequencies or higher-order super harmonics. 
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Table 1: CPU times for off-line and on-line computations. 



Numerical Experiment 

Procedure 

Airfoil 

Gavity, 

Low-Re 

Gavity, 
Moderate-Re 

FOM # of DOF 

360,000 

288,250 

243,750 

Time-integration of FOM 

7.8 hrs 

72 hrs 

179 hrs 

Basis construction (size n p ROM) 

0.16 hrs 

0.88 hrs 

3.44 hrs 

Galerkin projection (size n-\- p ROM) 

0.74 hrs 

5.44 hrs 

14.8 hrs 

Stabilization 

28 sec 

14 sec 

170 sec 

ROM # of DOF 

4 

4 

20 

Time-integration of ROM 

0.31 sec 

0.16 sec 

0.83 sec 

Online computational speed-up 

9.1 X 10^ 

1.6 X 10“ 

7.8 X 10“ 


Finally, a snapshot of the predicted velocity and pressure magnitudes at the final 
snapshot are illustrated in Figure 11 and Since the flow at this higher Reynolds 
number is chaotic, the low-dimensional model can not be expected to track the original 
snapshots exactly. However, the snapshots demonstrates that the stabilized ROMs faith¬ 
fully reproduce the large features of the flow. The same cannot be said of the standard 
ROMs. 


4-3. Computational speed-up 

For each problem considered, the speed-up factor delivered by its ROM for the online 
computations is reported in Table All ROMs are solved in MATLAB using the vari¬ 
able order solver 0DE15s; Jacobians are provided analytically. For more details on this 
algorithm, the reader is referred to Shampine et al. [43]. All online time-integration CPU 
times were measured using the tic-toc function on a single computational thread via 
the -singleCompThread start-up option. The FOM time-integration, basis construction 
and Galerkin projection times given in the table are reported in CPU-hours, calculated 
as the product of the number of processors used in the computation and the mean CPU 
time over all processors. The number of processors employed varied between 1 and 128. 
The online speed-up is calculated by evaluating the ratio between the time-integration of 
the FOM and the time-integration of the ROM. The reader can observe that the ROM 
online speedup is on the order of at least 10"^ for all three problems considered. Moreover, 
the stabilization step takes very little time (on the order of seconds/minutes). 


5. Conclusions 

In this paper, an approach for stabilizing and enhancing projection-based fluid ROMs 
of the compressible Navier-Stokes equations is developed. Unlike traditional approaches, 
no empirical turbulence modeling terms are required, and consistency between the ROM 
and the full order model from which the ROM is derived is maintained. Mathematically, 
the approach is formulated as a trace minimization on the Stiefel manifold. The method 
is shown to yield both stable and accurate low-dimensional models of several represen¬ 
tative compressible flow problems. In particular, the method is demonstrated on flows 


18 






















Figure 3: Domain and mesh for viscous channel driven cavity problem 


at higher Reynolds number where the dynamics are chaotic. Future work will include 
the extension of the proposed approach to problems with generic non-linearities, where 
the ROM involves some form of hyper-reduction (e.g., DEIM, gappy POD) following 
the procedure described in Remark 2, as well as to predictive applications with varying 
Reynolds number and geometry. 
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Figure 4: Nonlinear model reduction of channel drive cavity at Re 5^ 1500. Evolution of modal energy 
(a) and phase plot of the first and second temporal basis, ai{t) and a 2 (t) (b); DNS (thick gray line), 
standard n = 4 ROM (dashed blue line), stabilized n,p = 4 ROM (solid black line). Stabilizing rotation 
matrix, X (c) 



Normalized frequency (xtt rad/sample) 


Figure 5: PSD of p{x,t) where x = (2, —1) of channel drive cavity Re 5^ 1500. DNS (thick gray line), 
stabilized n,p = 4 ROM (black line) 
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Figure 6: Snapshot of channel drive cavity Re ^ 1500; contours of n-velocity magnitude at the final 
snapshot. DNS (top), standard n = 4 ROM (middle) and stabilized n,p = 4 ROM (bottom) 
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Figure 7: Nonlinear model reduction of channel drive cavity at Re 5^ 5500. Evolution of modal energy 
(a); DNS (thick gray line), standard n = 20 ROM (dashed blue line), stabilized n,p = 20 ROM (solid 
black line). Stabilizing rotation matrix, X (b) 
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Figure 8: PSD of p(xi, t) where xi = (2, —0.5) of channel driven cavity at Re Ri 5500. DNS (thick gray 
line), stabilized n,p = 20 ROM (black line) 



Figure 9: PSD of p{x 2 ,t) where X 2 = (0, —0.5) of channel driven cavity at Re Ri 5500. DNS (thick gray 
line), stabilized n,p = 20 ROM (black line) 
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Figure 10: CPSD of and p{x 2 ,t) where xi = (2,—0.5) and X 2 = (0, —0.5) of channel driven 

cavity at Re 5500. DNS (thick gray line), stabilized n,p = 20 ROM (black line) 
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Figure 11: Snapshot of channel drive cavity Re ps 5500; contours of w-velocity magnitude at the final 
snapshot. DNS (top), standard n = 20 ROM (middle), and stabilized n,p = 20 ROM (bottom) 
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Figure 12: Snapshot of channel driven cavity Re ps 5500; contours of pressure at the final snapshot. 
DNS (top), standard n = 20 ROM (middle), and stabilized = 20 ROM (bottom) 
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Appendix A. Construction of Galerkin matrices for the compressible Navier- 
Stokes equations 

In this section the construction of Galerkin matrices for the compressible Navier- 
Stokes equations Q is outlined. Consider the standard orthonormal POD vectors u*, 
n*, C*) and p*, * = 0, • • • ,n where * = 0 identifies the constant mean flow. The following 
products are generated for j, fc = 0, • • • , n: 


4^’^^ = Qul + e QpI- © 


.3 =u^Qv^,+v^(Dv^y+C^(Dp^y-—C^(D 

QpI+ + 7pJ' © {v^ + ul) - 


4 fe 2 

r" - r' 


- 3^^ 


+ {4 + 


+ {4 + 


(A.la) 


(A.lb) 
, (A.lc) 


7 


RePr 


[(p^ © Qk)xx ~\~ {p^ © Cfc) 


yy} 


+ 


1-7 

Re 


< © © oVy - 7 , 4 ] + (4 + © (4 + 


(A.ld) 


where subscripts identify partial spatial derivatives, and © is the Hadamard element-by- 
element product. A standard Galerkin projection for i, j,k = 0, - ■■ ,n is performed 
as follows: 


N 

{j, k)=Yl +u"Q + n* © -f p* © , (A.2) 


m—0 


where /i is a vector of element volumes. Finally, the standard Galerkin matrices for 
f, j, fc = 1, • • • ,n with oo = 1, are given by 


Q = 1>®(0,0), 

L,,, =6«(j,0) + 5W(0,j), 

Qfl=b^^Kj,k)- 


(A.3a) 

(A.3b) 

(A.3c) 
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